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ABSTRACT 

We modify an existing magnetohydrodynamics algorithm to make it more compatible with 
a dimensionally-split (DS) framework. It is based on the standard reconstruct- sol ve-average 

CN| strategy (using a Riemann solver), and relies on constrained transport to ensure that the mag- 

^^ netic field remains divergence-free (V • B = 0). The DS approach, combined with the use of 

C3 a single, cell-centred grid (for both the fluid quantities and the magnetic field), means that the 

^ algorithm can be easily added to existing DS hydrodynamics codes. This makes it particularly 

useful for mature astrophysical codes, which often model more complicated physical effects 
on top of an underlying DS hydrodynamics engine, and therefore cannot be restructured eas- 
ily. Several test problems have been included to demonstrate the accuracy of the algorithm, 

i— i and illustrative source code has been made freely available online. 

^h Key words: MHD - methods: numerical - magnetic fields 

• 

Oh 

q 1 INTRODUCTION 

When developing complex astrophysical codes, we often start with a simple hydrodynamics engine. This is designed to numerically solve 
^ the Euler equations of pure hydrodynamics, which - in conservative form - can be written as: 

dU dF dG dH_ 
£\\ dt dx dy dz 



^o 



% 



where U is the vector of state quantities we wish to evolve in time, and the flux vectors F(U), G(U) and H(U) contain their respective fluxes 



> 

On 

(-*) in the three spatial (Cartesian) dimensions. A popular approach for solving these equations numerically is the reconstruct- solve- average 
\^ (RSA) strategy (a term coined by Mignone et al. 2007). We start by dividing our space into a 3D grids of 'cells', and then we discretise the 
CO initial data by storing only the cell- volume averages of the state quantities at every cell centre. Each timestep then comprises three stages. 
<^|- First, we reconstruct this discretised data to estimate the value of the state quantities on either side of each cell face, thus defining - in 
C^) effect - a series of 'Riemann problems' r\ Next, we solve each of these problems independently using a Riemann solver algorithm, which 
CO determines the fluxes of the state quantities through each cell face. Finally, we use these inter-cell fluxes to determine the new cell-averaged 
\ I state quantities. These three stages are repeated until the desired amount of 'time' has been simulated. 

^ This strategy is often combined with dimensional splitting (DS), where the 3D grid is split into several independent ID problems and 

each timestep then becomes a series of 'directional passes'. In each of these passes, all of the rows in a given direction are evolved according 
to the RSA strategy above. However, we only deal with the cell faces which are perpendicular to the direction of the pass: the x-direction pass 
C^ determines the flux vector F(U) through all the x = constant faces, while the y- and z- direction passes determine G(U) and H(U) through 
the y = constant and z = constant faces respectively. The popularity of this approach is due to its many advantages: the timestep is not as 
constricted by stability conditions as with a fully multi-dimensional code; it is simple to upgrade an existing ID code to more dimensions in 
this way; the rows in each pass can be evolved independently, which means the code can be easily parallelised; it has been shown that the 
accuracy of the ID code can be preserved by alternating the order of the passes (Liska & Wendroff |2003| ). 

Once this hydrodynamics engine is in place, it can be improved by adding other physical phenomena (such as viscosity, self-gravity or 
a realistic equation of state), and more sophisticated computational features (such as parallelisation and adaptive/non-uniform grids). This is 
generally an effective way to develop code, but can cause difficulties if we then need to modify the underlying engine - for example, when 

* E-mail: hari.sriskantha@ed.ac.uk 

1 A Riemann problem is defined as the solution to a conservation law with piecewise constant initial data and a single discontinuity (LeVeque |l992) . 
Algorithms exist to find this solution both exactly and approximately (the latter sacrificing some accuracy in favour of computational speed). 
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(i*,j, k*) 
(i*J, k) 



Figure 1. This figure depicts the cell centred at position (i, j, k). The orange face is labelled (i*, j, k), where i* = i + 1/2, because it is located |Ai away 
from the cell centre. Using similar reasoning, the red edge is labelled (z*, j, k*). (See ^2.1| for details.) 

adding the effects of magnetic fields. The dynamics of a fluid in a magnetic field are governed by the equations of magnetohydrodynamics 
(MHD), which combine the Euler equations with Maxwell's equations of electromagnetics. Although these equations can be written in the 
same form as eq. ([TJ, they cannot be easily incorporated into DS-RSA codes because of the need to maintain the divergence-free condition: 

V • B = 0. (2) 

Ignoring this condition in numerical simulations can result in unwanted effects, such as the creation of unphysical forces parallel to the 
magnetic field, and the loss of conservation of momentum and energy (Balsara & Spicer [T999] ). 

There are several families of algorithms for maintaining this condition (which are summarised in e.g. |T6th|2000| ). In this paper, we focus 
on the constrained transport (CT) family, where the magnetic field quantities are discretised in such a way that V • B does not change in time 
analytically (to machine error). We then simply choose initial conditions such that V B = 0. Many CT algorithms have been developed (such 
as |Dai & Woodward|1998l|Balsara & Spicer|1999l|T6th|2000||Londrillo & del Zanna|2"000| ), but these are often dimensionally unsplit and/or 
require a 'staggered' grid: the fluid quantities are stored as volume- averages at cell centres, while the magnetic field quantities are stored as 
surface-avemges at cell faces. This is not troublesome for researchers intending to develop new codes from scratch, but is less than ideal for 
those working on mature codes as described above. Even implementing the non-staggered algorithm suggested in Toth (2000) - which was 
modified from the algorithm of Balsara & Spicer ( 1999} - would require considerable restructuring of the underlying hydrodynamics engine. 
Instead, a preferable algorithm is one that minimises the need for such restructuring. In this paper, we show that one can be developed by 
reframing Toth's approach so that it is more compatible with a dimensionally- split framework. By doing so, we aim to considerably simplify 
the development required to add the MHD equations to mature astrophysical codes. 

In §2 we describe some practical details, including notation, units and the formulation of the MHD equations used. In §3 we summarise 
the Balsara-Spicer-Toth algorithm for maintaining the divergence-free condition, while in §4 we derive our 'more' dimensionally- split ver- 
sion. Readers who are only interested in the final algorithm may skip ahead to §5, where we include a clear step-by-step guide to implementing 
it. (Additionally, illustrative source code has been made available online - which was designed to be as easy as possible to be read - in order 
to aid understanding of the algorithm.) Finally, in §6 we demonstrate its accuracy with several well-known MHD tests. 



2 FORMULATION AND NOTATION 

2.1 Discretisation 

We model Cartesian space in the domain D x x D y x D z . We divide this space into n x x n y x n z 'cells' in the usual way, each of volume 
Ax x Ay x Az, such that D x — n x • Ax, etc. (We assume that the cell sizes are constant and uniform, but it is simple to extend our 
algorithm to adaptive/non-uniform grids.) The position of each cell centre is labelled as (i, j, k), where i E {1, 2, ..., n x }, etc. We also refer 
to the faces and edges of cells, using starred notation i* = i + 1/2. A face will therefore have one starred co-ordinate, e.g. (i*,j, k), while 
an edge will have two, e.g. (i*,j*,h). A typical cell, face and edge are depicted in Fig. [I] Finally, we discretise time such that t n = n- At n . 
The timestep At n can change with time, but must always satisfy the Courant-Friedrich-Lewy (CFL) stability condition. 

To avoid excessive indices, we only write down spatial co-ordinates that differ from (i, j, k). For example, the face-average Un ,j,k+1 
will be abbrievated to t/^*' fe+1 . The reader should therefore assume that there are always three spatial co-ordinates, even if all three are not 
explicitly written down. When there is no change from (i,j, k) we use *, such that U\i j,k is abbrievated to £/*. Similarly, we do not use 
subscripts to denote components of discretised values. For example: the cell average of the x-component of the magnetic field will be written 
as (Bx)n, not B xn . A summary of the non-standard notation used throughout this paper is provided in TablefT] 

Using the discretisation above, we can now write a numerical equation for the DS-RSA strategy, which gives t/* +1 in terms of U* n : 

77* — 77* - —(F 1 * - F i *~ 1 \ - —fr j * - r j *~ X \ _ ^L(fj k * _ l^* -1 ^ n\ 

V n-\-l V n i, f n* r n* ) A V^^n* ^n* ) A \"-n* "n* )i \ J ) 

Zax l\y l\z 
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Table 1. A summary of the non-standard notation used throughout this paper. See §2.1| for more details. 



Notation 



Meaning 



the cell centred at (i, j, k) 

the cell centred at (i, j + 1, k) 

the index i + 1/2 

the cell face centred at (i*, j, k + 1) 

inter-cell flux vectors returned by the Riemann solver 

the 7th component of F 

the x-component of B 

cp = 1 or 1/a/4tt depending on units (see { 2.2 i 



a + 1) 

(z*,fc + 1) 
F,G,H 

(Ft) 
(Bx) 



where the last three terms correspond respectively to the three directional passes. The terms F nsle , G n * and H n * are the inter-cell flux vectors 
returned by the Riemann solver during the 'solve' stage (signified by tildes), and n* = n + 1/2 refers to the half-timestep. 



2.2 Units 

In the literature, there are two sets of units commonly used for the magnetic field, which differ by a factor of \^4tt. To aid compatibility, we 
use unit-neutral notation, using a constant ip. For SI units, ip = 1, while for CGS units, ip — l/y^r. By including this constant in the code, 
we can switch easily between the two. We will make it clear which version is used by any referenced papers. 



2.3 The MHD equations 

The MHD equations can be written in the same form as eq. ([TJ. The vector of state quantities is defined as: 
U= [p,pv x ,pv y ,pv z ,£,B x ,By,B z ] T ', 



(4) 



where p is the mass density, pv is the momentum density, S is the total energy density, and B is the magnetic field density. Recall that the 
cell-averaged values of these quantities are stored at every cell centre. We refer to the first five elements as 'fluid quantities' and the last three 
elements as 'magnetic field quantities' whenever we need to distinguish between the two. The flux vectors are defined as: 



F = 



pv x 

2 i 2 7->2 

pv x + p T - <p B x 

PV X Vy - (f 2 B x By 

pv x v z - ip 2 B x B z 

v x (£+p T )-(p 2 B x (vB) 


-E z 
E v 



G = 



pVy 
pV X Uy - (f 2 B x By 

PVy +PT -¥ 2 By 

pVyV Z - if 2 ByB Z 

Vy(S+ PT )-(p 2 By(vB) 

E z 


—E x 



H = 



pv z 



pv x v z 

PVyV z 



ip 2 B x B z 

(f 2 ByB Z 



pv 2 +pt - p 2 B 2 
v z (£ + p T )-(p 2 B z (vB) 

~Ey 

E x 




(5) 



where px is the total pressure, such that px = Pm + Pg- The magnetic pressure is given by pm = \tp 2 B 2 ' , where B 2 
pressure, pc is related to the total energy density by an ideal equation of state, which closes the system, such that: 

c PG . 1 2 . 1 2 D 2 

where 7 is the ratio of specific heats. Finally, E is the electric field, which is given by the ideal version of Ohm's law: 
E= -v xB. 



B • B. The gas 



(6) 



(V) 



3 BALSARA-SPICER-TOTH ALGORITHM 



We now describe the algorithm of Balsara & Spicer ( 1999), later modified by Toth (2000). The fluid quantities are evolved using the usual 
DS-RSA strategy, as described in §1, while the magnetic field quantities are evolved using constrained transport. 



3.1 Constrained transport (CT) 

CT was first developed by Evans & Hawley (1988), and uses Faraday's law to determine how the magnetic field B changes in time: 

dt B 



-V xE. 



(8) 
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where E is the electric field. If take the divergence of both sides of this equation, then we get: 
d . 



dt 



V B = -V- (V xE) = 0. 



(9) 



where the second equality is true because V • ( V x F) — for any vector F. This means that if we evolve the magnetic field using eq. {8f , 
then V • B will not change in time analytically. We can therefore maintain the divergence-free condition by simply choosing initial conditions 
such that V • B — (Toth 2000). To implement eq. (8) numerically, we start by considering the x-component of the magnetic field averaged 
on the face (i*), which we call (foe)**. If we define (i*) as having area AA, then it follows that: 



d_ 

dt 



(bxy 



AA 



V xE-dA. 



(10) 



Note the use of the lower case b to denote a face-averaged magnetic field quantity. (We will later use the upper-case B to distinguish the 
ce//-averaged magnetic field quantities.) Then, using Stokes' theorem, we have: 



d ( w* 1 

— (ox) = — — — - 
dt v ; AA 



E • dr. 



(11) 



where we now integrate over the path dA, which is the boundary that surrounds the face (i*). We can discretise this equation, because we 
know this boundary comprises four edges in discretised space: (i*, j*), (i*, j * —1), (i*, /c*) and (i*, k * —1). Furthermore, if we integrate 
in time, and use the fact that AA = AyAz, then we obtain the following expression for (frx)^* +1 given (6x)^*: 



At 



At 



[(E Z y:r - [Ez^r- 1 ] + ^ [(E y y: : k * - (Eyy:^- 1 ] , 



\OX ) n -\-l \OX) n ^y-^f^jn* y-^,~j n * j , 

where (Ez) l * ,j * is the z-component of the electric field, but averaged on the edge (i*,j*): 



(Ez. 



i*,j* __ I / 

Az , . 



E-dz. 



(12) 



(13) 



(i*,j*) 



By substituting eq. |T3| ) into eq. |T2l >, it should be clear how the latter follows from eq. (\J\ . Using similar reasoning, we can then derive 
expressions for the other two (face- averaged) components of the magnetic field, (6y)^* +1 and (bz)n* +1 : 



(by)^ = (Wn 



M [(ExYZ^ - (Ex)^' 1 ] + ^ [{Eztr - (Ezft- 1 **] 



(bz)n+l 



(bz) ri 



Az 
At 
Ax 



[(Ey)l 



{Ey)l 



Ax 
At 



] + ^ [(Ex)^* - (Ex)ir 1M ] • 



Ay 



(14) 
(15) 



3.2 The electric field 



The next step is to determine the electric field components in eqs. ( 12 1, ( 14 >-< 15 L Balsara & Spicer |l999 (p — l/V^r) use the fact that 
these components appear in the last three rows of the MHD flux vectors - given in eq. {5} - as fluxes of the magnetic field quantities. This 
means that face-averaged values of these components will be returned by the Riemann solver during the solve stage. We can convert these 
into the edge- averaged values we require by simple arithmetic averaging; to make the explanation of this clearer, we re-label the last three 
rows of the flux vectors as follows (where tildes signify that these are the values returned by the Riemann solver): 



F : 








' (G6)& " 




\ (H6)t 


(F7)Z 


, G: 





, H. 


(H7)t 


(Fsy n x 




(Gsyz 








(16) 



So for example, (Ex) appears as (GS) 3 n* on y = constant faces and as (H7)nl on z = constant faces. In our grid of cells, we know that each 
edge is adjacent to four faces. For example, the edge (j*, k*) is adjacent to the faces (k*), (j + 1, k*), (j*) and (j*, k + 1), as depicted in 
Fig.[2p. So to obtain the x-component of the electric field averaged along this edge, (Ex) 3 n* *, Balsara and Spicer simply take the arithmetic 
mean of the appropriate fluxes through the surrounding four faces (where minus signs compensate for the minus signs in eq. J5])): 

(£*)&*• = \ [{H7) k n l + {H7Y+ 1M - (G8)H - (G8)i'- k+1 ] , (17) 

and similarly for the other components: 

(EyW = - [(F8)t + (F8y:: k+1 - (h 6)*: - (my+ 1M ] , 



{Ezy:: k * = \ [(G6)t + (G6);t lj * - (F7y n i - (F7y n y +1 ] 

This completes the description of the Balsara-Spicer algorithm. 



(18) 



(19) 
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(i* k*) 



(j* k*) 



/ 



U + 1. k*) 



V 

U*) 



(a) 




(i* k*-l) 



(i*,j*-l) 




(i*,j*) 



(b) 



Figure 2. (a) The edge-averaged electric field. The edge (j*, k*) is highlighted in red. To determine the edge average of (Ex) along this edge, we take the 
arithmetic mean of the face averages of (Ex) on the surrounding four faces. These, in turn, are returned by the Riemann solver as inter-cell fluxes during the 
solve stage. See §3.2| for more details, (b) The face-averaged electric field. We can then determine two different face-averaged components of the electric field 
on each face. For example, on the face (z*) - highlighted in green - we can determine (Ey) l _^ x by avera ging (Ey) along the edges (z*, k*) and (z*, k * — 1) 
(top left), or (Ez)% x by averaging (Ez) along the edges (z*, j*) and (z*, j * — 1) (bottom right). See j|3.3|for more details. 



3.3 The cell-averaged magnetic field 

It is clear that the Balsara-Spicer algorithm requires a staggered grid: the magnetic field quantities are evolved on the cell faces, while the 
fluid quantities are evolved at the cell centres. Whenever a cell-averaged magnetic field is required (for example, during the reconstruction 
stage), then we must interpolate. A linear interpolation is usually adequate: 



(Bx) % v 



- [(bx)n + (bx)n *] 



(20) 



where we now use the upper-case B to denote a cell-averaged magnetic field. Toth (2000) proposed a simple extension to this idea: to build 
the interpolation into the algorithm itself, thus eliminating the need for a staggered grid. His paper only includes the equations for a 2D 
system, but as he mentions, it is straightforward to extend these to three dimensions: 



{Bx)' n+1 = (Bx)' n - -g [{Ez)% - [EzY*- 1 ] + ^ [(Ey)*, - {Ey)^ 1 ] 



{By)l- 



(Bz)' n 



n+1 



(By)l 
(Bz)' n 



At 

Az 
At 
Ax 



[{Ex) k . 
[(Eyt 



{Ex)^: 1 ] + 



(Ey)- 



Az 

At 

Ax 
At 



[{Ezf. 



+ ^ [^Y- 



(Ez) 



(Exy- 1 ] 



(21) 
(22) 
(23) 



which is similar to eqs. (T2) , fj"4)-fi"5}, except we now require /ace-averaged electric field components. Note that these are not the same 
as the face-averaged components we obtained earlier from the Riemann solver. Instead, they are found by taking the arithmetic mean of 
the edge-averaged electric field components - of eqs. fl7)-([T9} - in two surrounding, parallel edges. Each face has two such face-averaged 
electric field components, which depend on the two edges we use to calculate this mean. For example, for the face (i*), we have: 

{Ey)% x = \ [(Eyt> k * + {Ey)"' 1 "- 1 ] and {Ez)% x = \ [{Ezf^ + {Ez) 1 *' 3 '- 1 ] . (24) 

Both of these are illustrated in Fig.J2p. We use the — >> subscript to distinguish, for example, between (EyY_X x through the x = constant faces 
and and (Ey)^ z through the z = constant faces. Using eqs. {l7|l-( 19 1, we then have: 



{e v )% x = x - [{Fsy:; k + {F8y:; h+1 - {my** - {my+ hk * + {Psy:*- 1 + {Psy:* - {my**- 1 - {my+ hk *- 1 ] , (25) 

{e z )% x = | [{G6y n r + (G6);t lj * - {P7y n y - {P7y n y +1 + {G&U*- 1 + {G^y*- 1 - {P7y n y~ 1 - {P7y n r] , (26) 

In a similar manner, we can derive the face-averaged components through the faces (z*) and (&*): 

{Ex)% = i [{Hiy** + {H7)it 1M - {G8)i*' k - {G8)£' k+1 + {my**- 1 + (HT) i + 1 ' k *- 1 - (GS)^" 1 - {G8)i*' k ] , (27) 

{ezy: v = I [(G6)Ui* + {Gey n + y* - {P7y n y - {P7y n y +1 + {G6y- hJ * + {gqW - {P7y n *~ h3 - (ft)*- 1 ** 1 ] . (28) 

[Ex) k 4 x = \ [{H7) 3 ** + {H7)it 1M - {G8)i** - {G8)i** +1 + {HTfc 1 - km + {H7) 3 ** - {G8)iT 1 ' k - (G8)ir M+1 ] , (29) 
© 2013 RAS, MNRAS 000,[T]{T2] 
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(Ey) k _: z = | [(P8y:: k + (Fsy:: k+1 - (my** - (my+ 1M + (FsjJrr 1 -* + (F8y n r hk+1 - (my n ~ hk * - (my**] . po) 

This concludes the description of the Balsara-Spicer-Toth algorithm. 

3.4 Disadvantages of this algorithm 

The simplest implementation of this algorithm involves creating six new 3D arrays, one for each of the inter-cell fluxes we require: (F7), 
(F8), ((56), ((58), (H6) and (HI). Every time we complete the solve stage for a given ID row, we obtain the values of two of these fluxes 
along the entire row, and so can store them in the appropriate arrays. (The two we obtain depend on what the directional pass is.) Once a 
complete set of three directional passes has been completed, we can sort these arrays into another six 3D arrays, one for each of the face- 
averaged electric field components: (Ex) J J^ y , (Ex)^ z , (Ey)™ x , (Ey)^ z , (Ez) z ^ x and (Ez) j ^ y . This is done according to eqs. (25l-(30l 
above. Finally, we use these components to update the magnetic field quantities using eqs. (21]) -(23]). 

This algorithm has been shown to be reliable through standard MHD tests (Toth 2000 ), but has some notable disadvantages: 

(i) Storage of information. In a DS-RSA hydrodynamics code, we can extract a single row from the 3D grid, update all of the quantities 
along that row, and then return it. We do not need to store any information from that row to evolve other rows in the same directional pass. 
This is clearly not the case here: we need to store all of the inter-cell fluxes and electric field components. 

(ii) Processing between sets of directional passes. In a DS-RSA hydrodynamics code, once a complete set of directional passes has been 
completed, the state variables require no further operations until the next set. Again, this is not the case here: we need to carry out additional 
operations on the inter-cell fluxes and electric field components to determine the new magnetic field quantities. 

(iii) Staggered grids. Although the magnetic field is stored as cell averages, we still require staggered grids: the inter-cell fluxes and 
electric field components are stored as face averages. This complicates their implementation, especially with adaptive/non-uniform grids. 

(iv) Limited parallelisation. Finally, because of the extra processing between directional passes, and the fact that the electric field compo- 
nents depend on inter-cell fluxes from different rows and directional passes, it is not immediately obvious how this code could be parallelised. 

While these are arguably not major issues, they complicate the implementation of this algorithm in existing, mature astrophysical codes. 
We would therefore prefer an algorithm that: minimises the amount of data that needs to be stored from each row during the directional passes; 
reduces the amount of processing required between complete sets of directional passes; completely removes the need for any staggered grids; 
simplifies and maximises the possible parallelisation. We now show that such an algorithm exists, by reformulating the Balsara-Spicer-Toth 
algorithm so that it is more compatible with a dimensionally- split framework. 

4 A MODIFIED APPROACH 

To start, we define just three new 3D arrays: A(Bx), A(By) and A(Bz). At the end of each complete set of three directional passes, we aim 
for these to contain the change in the magnetic field during those passes. To make our derivation clearer, we focus on a single cell: (i, j, k). 

4.1 The x-direction pass 

4.1.1 Contribution to A(Bx) m 



The change in A(Bx) m during each timestep is given by eq. (21 1, and depends on the components (Ez)iX y , (Fzyj^y 1 , (Ey)^ z and 



(Ey)^ z 1 . We can substitute in the definition of these components, given by eqs. (30 1 and (28 1. However, we only include the F terms, as 



these are the only ones we obtain during the x-direction pass. Then, after simplification, we are left with the following expression: 
A(Bx)' := A(Bxy +£^§ [(F7)^' +1 + (F7)£r 1J+1 - (i^)^" 1 - (F7)t" 1J_1 ] 

+^k[(F8y:i k+l + (F8y:r 1 ' k+1 -(F8y:i k - l -(F8y n r 1 ' k - 1 ], ( ' 

where use ':=' to refer to an update in the usual computational sense. 

This is still not in the spirit of dimensionally- split algorithms, as the right-hand side of the expression depends on fluxes from several 
different rows, as shown by the different spatial indices. Instead, we want to consider what we could do after a given x-row - labelled 
(i, J, K) with fixed J,K - has been evolved. To do this, we use the fact that if A( Bx) m depends on, say, (F7)n*' J+1 , then this is equivalent 
to saying that A(Bx) j ~ 1 depends on (F7)n*. If we consider each term in eq. (3l) in this manner, then we can rewrite it as: 

A(Bx) J±1 := A(Bx) J±1 T ^ | [{FT)* + (FT)"- 1 ] , (32) 

A(Bx) K±1 := A(Bx) K±1 T ~ [(F8)£ + (FS)^ -1 ] . (33) 

Both the right-hand sides only contain terms from the row we are currently on, (i*). This means that once a given row has been evolved, we 
can process the contribution of the inter-cell fluxes to the magnetic field quantities immediately - using eqs. (32) and (33) - without needing 
to store them in any other form. It is also worth noting that these equations are near identical to the standard RSA approach given in eq. (3). 

© 2013 RAS, MNRAS 000,[l}[T2] 
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Table 2. In a dimensionally -split code, we can avoid rewriting the ID algorithms if we make sure that the Cartesian components of the velocity are extracted 
(and then returned) according to this table. A similar strategy applies for the components of the magnetic field. 



Dimensional pass x 



(vx) 1B stores (vx) (vy) (vz) 

(t;y) 1D stores (vy) (vx) (vx) 

(vz) 1B stores (vz) (vz) (vy) 



4.1.2 Contribution to A(By) m and A(Bz) m 

Following a similar approach, we use eq. (J22| to find the contribution to A(By) m during the x-direction pass: 



A(fly)' := A(fly)' + ~ [-(F7)* + (FT)^ 1 ] , (34) 

AiBy)^ 1 := A(By) J±1 + ~ [-(F7)£ + (F^f 1 ] ■ (35) 

while eq. ( [23] ) gives us the contribution to A(Bz) m during the x-direction pass: 

A(BzY := A(ifc) - + £* i [-(F8)£ + (FS)^ 1 ] , (36) 



Ax 4 
z 
Ax; 



A(^) fc±1 := A(Bz) k±1 + ~ [~(F8)Z + (FS)ST 1 ] • (37) 



4.2 The 2/- and z-direction passes 

We go through the same steps as above for the ^/-direction pass, but start by considering A(By) m . The change in this during each timestep 



is given by eq. (22i, and depends on the components (Ex) k X z , (Ex)^ z 1 , (Ez)% x and (Ez)^^ . Again, we substitute in the definition of 
these components, given by eqs. J29| ) and J26| , although this time we only include the G terms. After simplification, we are left with: 

A(fly)' := A(By)' + £| [{GQ)^* + (GO)^*" 1 - (G6)^* - (Gefc 1 '''- 1 ] 

+|f [(G8)^' fc+1 + (G8)ir 1,fc+1 - (GSYnV"- 1 - (GS)^" 1 '*- 1 ] . l J 

Before we continue, a note on dimensionally-split codes: when we extract a given row from each full 3D array (of which there is one for each 
state quantity), we store them in corresponding ID arrays. However, the quantities that are broken down into Cartesian components - the 
velocity and the magnetic field - must be extracted differently depending on the directional pass. For example, in the x-direction pass, (vx) is 
parallel to the row, and is stored in a ID array we shall call (vx) 1B . However, in the ^/-direction pass, it is (vy) that is parallel to the row. To 
avoid having to rewrite the ID algorithms for each directional pass, we simply store (vy) in the array (vx) 1B . This approach is summarised in 
Table[2] If we rewrite the terms in eq. J38l according to this approach - such that {x, F, i, 6} — >• {y, G, j, 7} and {y, G, j, 7} — ► {x, F, i, 6} 
- then we actually retrieve eq. ( [37] ). Furthermore, it can be shown that all of eqs. ([32|-([37]) can be retrieved in this manner, and for the 
z-direction pass as well. In other words, once the equations for the x-direction pass has been implemented, we do not need to re-write them 
for the y- and ^-direction passes, as long as the Cartesian components are extracted and returned as described above. 

4.3 Final update 

Once all three passes have been completed, we simply need to perform the following operation: 

(Bx) := (Bx) + A(Bx) (39) 

and similarly for (By) and (Bz) to update the magnetic field. There are some additional subleties associated with this, regarding the pressure 
and energy (which depend on the magnetic field), and these are described in the next section. 

4.4 Comparison with the Balsara-Spicer-Toth Algorithm 

(i) Reduced storage of information. We have minimised the amount of information that needs to be stored during the directional passes. 
Instead of storing the individual inter-cell fluxes and electric field components (which requires twelve 3D arrays), we only store the change 
in the magnetic field quantities (which requires just three). 

(ii) Reduced processing between sets of directional passes. We have also minimised the number of operations required between complete 
sets of directional passes: we only need to perform the simple calculation given in eq. ( [39] ). 

(iii) No staggered grids. Furthermore, the three arrays we require store only c^//-averaged values, which means we no longer need to store 
any information on the cell faces. This simplifies the use of adaptive and non-uniform grids. 
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(iv) Simplified parallelisation. The modified algorithm is more easily parallelisable. For example, assume that we have a number of 
processors, each independently evolving some rows in a given directional pass. Each of these processors should have a copy of the AB array, 
which are updated as described above. Then, because these all contain changes in the magnetic field quantities, we only require a simple 
'collect' -style operation to sum over all of these copies and obtain the total change before completing eq. (39}. 

(v) Dimensionally -split framework. Finally, because this algorithm is more compatible with the dimensionally- split framework, it is con- 
siderably simpler to include in existing dimensionally-split hydrodynamics codes. 



5 THE COMPLETE ALGORITHM 

We will now describe the complete algorithm for each timestep. Note that we require the following sets of arrays: 

(1) 7 x t/ 3D : the 3D arrays which store the cell-averaged state quantities as listed in eq. ^. 

(2) 7 x £/ 1D : the ID arrays which are used to store ID rows extracted from the full 3D grid. 

Additionally, we require two more sets of arrays for the magnetic field quantities: 

(3) 3 x B s : used to store a copy of the cell-averaged magnetic field quantities at the start of each timestep. 

(4) 3 x AB: used to store the change in the cell-averaged magnetic field quantities during each timestep. 



5.1 Before the directional passes 

(a) Make a copy of the magnetic J 
corresponding B s arrays. These copies should not be modified during the directional passes. 



(a) Make a copy of the magnetic field quantities. Copy the three t/ 3D arrays that store the cell-averaged magnetic field quantities into the 



5.2 The x-direction pass 

(b) Extract a ID row. Choose an x-direction row from the grid, which is defined as having fixed j and k indices. Extract the cell-averaged 
state quantities along this row from the £/ 3D arrays, and store them in the corresponding t/ 1D arrays. 

(c) Reconstruct the cell-averaged data. For each U 1B array, interpolate the cell-averaged data to obtain left and right states at each cell 
interface (z*). These will be input states for Riemann problems. For the tests in this paper, we use the TVD scheme of Balsara ( 199 8] (p = 



Colella & Woodward 



lv47r), which includes a steepening algorithm. However, a popular alternative is the Piecewise Parabolic Method of 
{1984}, an MHD version of which is summarised in §4.2.3 of |Stone et aT] ( [2008| <p = 1). 

(d) Find the inter-cell fluxes. Solve the Riemann problem at each cell interface, to obtain the inter-cell flux for each state variable. We use 
the multi-state HLLD approximate Riemann solver of Miyoshi & Kusano ( 2005[ (p — l)[j 



(e) Update the state variables. Using the inter-cell fluxes F which were determined in step (d), update the cell averages stored in t/ 1D 



using eq. |3}. Note that this should be done for the magnetic field quantities as welljj 

(f) Return the ID row. Return the new cell averages stored in the t/ 1D arrays to the I 

(g) Complete the x-directionpass: repeat steps (b) to (f) for each ID row in the x-direction (that is, for all fixed j, k). 



(f) Return the ID row. Return the new cell averages stored in the t/ 1D arrays to the f/ 3D arrays (again, including magnetic field quantities). 



5.3 The y- and ^-directional passes 

(h) Complete all passes: repeat the equivalents of steps (b) to (g) for the y- and z- directions. Remember to extract and return the velocity 
and magnetic field components in the usual dimensionally-split manner, as described in Table|2] Also note that in some code implementations, 
the order of the passes is alternated between timesteps to preserve the accuracy of the ID method. 



2 A note on using non-staggered grids: in a ID system, the parallel component of the magnetic field - which is (Bx) in the ^-direction pass - should be 
uniform and constant in time. This automatically satisfies the divergence-free condition. The reconstruction and solver algorithms we use were designed for 
such ID systems, which means that the former does not find left and right states for this parallel component in step (c), while the latter assumes it is constant 
across the cell interface in step (d). This is not an issue when using staggered grids, because this constant value is stored explicitly at the cell interface. However, 
this is clearly not the case for non-staggered grids, which means we must interpolate this parallel component to get its value across the interface. Our tests 
have shown that a simple linear interpolation is sufficient, for example: (-B#)i n t,i* — [(Bx)i + (Bx)i+i]/2. 

3 This is because we require the gas pressure pq when we return to the earlier steps. However, recall that we only store the total pressure pt, so we need to 
calculate it by subtracting the magnetic pressure pm (i-SvPG — Pt — PmX which is in turn determined by the magnetic field quantities. Hence, in cases where 
Pm ^> PG> n °t updating the magnetic field components can result in negative pressures. This issue is discussed further in Balsara & Spicer||1999|. 
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5.4 After the dimensional passes 



After all of the above steps have been completed, the fluid quantities will be fully evolved to the next timestep. We will also have three 
versions of the magnetic field quantities: their values at the start of the timestep (stored in IT), their (correct) change during the timestep 
as determined by the constrained transport algorithm (stored in AB), and their values as evolved - incorrectly - through the standard RSA 
strategy (stored in the t/ 3D arrays corresponding to the magnetic field quantities). Then the final steps are: 

(i) Correct the energy and pressure (I). This is an optional step. We use the magnetic field quantities stored in the t/ 3D arrays to calculate 
the magnetic field pressure pm (using pm = \tp 2 B 2 ). We then subtract this contribution from the newly evolved values of the total energy 
(using the equation of state defined in eq. {6} as a guide) and the total pressure. 

(j) Update the magnetic field. We replace the magnetic field quantities stored in the U 3r> arrays with B s + AB (in other words, their 
value at the start of the timestep, plus their correct change in value during the timestep, according to the constrained transport algorithm). 
Consequently, the stored magnetic field quantities now maintain the divergence-free condition. 

(k) Correct the energy and pressure (II). If step (i) was completed: use these new magnetic field quantities to recalculate the magnetic 
pressure and add this back to the total energy and total pressure. This correction means that energy may not be conserved entirely (as 
discussed in Balsara & Spicer 1999), but avoids negative pressures if pm ^> Pg, as described in footnote[3] 

6 MHD TESTS 

To illustrate the accuracy of the modified algorithm, we now present the results of four standard MHD tests. All results were obtained using 
the second-order TVD reconstruction of Balsara (1998), and the HLLD Riemann solver of Miyoshi & Kusano (2005 ). Also note that we use 
units of if = 1/V^7r, except for the circularly -polarised Alfven waves, which use ip = 1. 

6.1 Circularly-polarised Alfven waves 

The test of circularly polarised Alfven waves were first described by Toth (2000). Such waves are smooth, analytic, non-linear solutions of 
the MHD equations, which makes them ideal for testing the implementation of these equations. The version implemented is sinusoidal and 
uses a periodic domain that can fit exactly one wavelength. We test two versions: a travelling wave, where the wave returns to its original 
position at integer multiples of the period (which is equal to 1 in this case), and a standing wave, where it stays at its original position. 

We use the initial conditions described in Toth (2000). We start by defining a rotated coordinate system, cy and c±. These are rotated 



by an angle = tt/6 from the Cartesian x and y (and hence our discretised grid) such that: 

C|| = xcosO + ysinO and c± = — xsinO + ycosO. (40) 

We then define the velocity and magnetic field components along these axes: v\\ = for travelling waves and v\\ = 1 for standing waves; 

B\\ — 1; v± = B± = 0.1sin(27TC||). We then have: 

(a) Space and time. Domain: < x < 1/ cos 0, < y < sin 0; number of cells: varying from 8 2 to 32 2 ; total time t — 5. 

(b) Fluid quantities. Across the entire grid: mass density p = 1; gas pressure pc = 0.1; x-velocity v x = v\\ cos — v± sin 0; ^/-velocity 
v y = v\\ sinO + v± cosO; z-velocity v z = 0.1 cos(27TC||); adiabatic index 7 = 5/3. 

(c) Magnetic field quantities. The magnetic field is set up equivalently to the velocity. Across the entire grid: cpB x = B\\ cos — B± sin 0; 
cpBy = B\\ sinO + B± cos 0; cpB z = 0.1 cos(27TC||). 

The results for both the travelling and standing waves at time t = 5 are shown in Fig. [3] and can be compared with Figs. 8 and 9 in 
Toth (2000). Note that we plot B±, which can be determined from B x and B y using eq. ( [40} . We include reference plots, which are simply 
the initial conditions on a grid with 128 x 128 cells. This was designed as a difficult test - especially at low resolutions - so the results do not 
coincide exactly with the expected analytic solution. However, they are comparable with results from other algorithms, as presented in |T6th| 
(2000). It may also be possible to improve them with the use of a higher-order reconstruction. 

6.2 The blast problem 

The blast problem is initialised as a small disk of high-pressure fluid in the centre of the grid, which then expands rapidly into the low-pressure 
ambient fluid. As pointed out by Balsara & Spicer (1999), it is not necessarily useful for testing whether the magnetic field is divergence 
free: the conditions are such that the build-up of non-zero V • B will not have a major effect on the overall dynamics. However, it is useful 
for testing the propagation of strong MHD shocks in multiple dimensions ( [Stone et al.|2008| ). 

We use the initial conditions described in Stone et al. (2008 ) (which were in turn taken from Londrillo & del Zanna 2000): 

(a) Space and time. Domain: -0.5 < x < 0.5; -0.75 < y < 0.75; number of cells: 300 x 450; total time t = 0.2. 

(b) Fluid quantities. Across the entire grid: mass density p = 1; velocity v = 0; adiabatic index 7 = 5/3. The gas pressure is set to 
pc = 10 within a radius r < 0.1 in the centre of the grid, andpc = 0.1 everywhere else. 
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(a) Results for travelling waves, v\\ = 0. 
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(b) Results for standing waves, v\\ = 1. 



Figure 3. Plot of B± = (V3B y — B x )/2, from the circularly-polarised Alfven waves test at time t = 5, initialised as described in 1 6.1 




(a) Mass density. (b) Thermal pressure. (c) Magnetic pressure. 

Figure 4. Results from the blast problem at time t = 0.2, initialised as described in j|6.2| 



(c) Magnetic field quantities. Across the entire grid: cpB x = cpB y — l/\/2; (fB z — 0. 

The mass density, thermal pressure and magnetic pressure at time t — 0.2 are shown in Fig. [4] and can be compared with Fig. 28 in 
|Stone et aL] ( [20081 ). The results are as expected: the blast - which would otherwise be circular - has aligned itself with the magnetic field. The 
near-perfect 7r-rotational symmetry of the results are also of note. 

6.3 The Orszag-Tang vortex 

The Orszag-Tang vortex is another commonly-used 2D MHD test, which starts with smooth initial data but becomes progressively more 
complex (Toth 2000 ). The expected symmetry of the solution is useful for showing that the dimensionally- split algorithm (which uses the 
same code for each dimensional pass, but just extracts and returns the rows differently) works and has been implemented correctly. 
We use the initial conditions described in Stone et al. (2008 ) (which were in turn taken from Orszag & Tang 1979| >. 

(a) Space and time. Domain: < x, y < 1; number of cells: 128 2 ; total time t = 0.5. 

(b) Fluid quantities. Across the entire grid: mass density p = 25/367r; pressure p = 5/12tt; x-velocity v x — — svn(2^y)\ ^/-velocity 
v y = sin(27rx); z-velocity v z = 0; adiabatic index 7 = 5/3. 

(c) Magnetic field quantities. Across the entire grid: cpB x — — sin(27ry)/\/47r; cpB y = sin(47rx)/\/47r; B z — 0. 



The mass density, thermal pressure and magnetic pressure at time t = 0.5 are shown in Fig. [5] and can be compared with Fig. 22 
in |Stone et al.| p008). The mass density can also be compared with Fig. 16 in |T6th| f2000| ) or Fig. 14 in |Dai & W oodward (1998J). The 
7r-rotational symmetry of the results implies that our dimensionally- split approach is effective. 
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(a) Mass density. (b) Thermal pressure. (c) Magnetic pressure. 

Figure 5. Results from the Orszag-Tang vortex at time t = 0.5, initialised as described in §6.3| 
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Figure 6. Results from the rotor problem at time t = 0.15, initialised as described in §6.4| All contours are linearly-spaced from 0.426 (blue) to 12.4 (red) for 
the mass density, 0.0369 (blue) to 1.98 (red) for the thermal pressure and 0.0154 (blue) to 2.60 (red) for the magnetic pressure. 



6.4 The rotor problem 



Finally, we consider the rotor problem, first proposed in its current form by Balsara & Spicer (1999). This involves placing a rotating disk 
in the centre of a square domain, which creates rotational discontinuities and propagates Alfven waves into the ambient fluid. The magnetic 
field is initialised such that it wraps around the rotor, squeezing the fluid from its original circular form. This problem is very useful for 
testing whether the magnetic field is divergence free, because the build-up of non-zero V • B = greatly influences the dynamics. 

We use the initial conditions described in |Balsara & Spicer| ( |1999| >. However, as noted by Toth (2000), the results they present are 
(erroneously) from a different set of conditions. We stick to what is described in the text, and compare our results with the correct results 
presented in other papers. We define the radius of the rotor in the centre of the grid as r < ro = 0. 1. This is separated from the ambient fluid 
by a taper, ro < r < n = 0.115. We also define vo = 2, and a constant / = (n — r)/(ri — ro) (which is used for the taper). Then: 

(a) Space and time. Domain: —0.5 < x, y < 0.5; number of cells: 256 2 ; total time t = 0.15. 

(b) Fluid quantities. The ambient fluid: mass density p = 1; pressure p = 1; velocity v = 0; adiabatic index 7 = 1.4. The rotor is defined 
within a radius r < 0.1: mass density p = 10; x-velocity v x = —vo(y — 0.5) /ro; ^/-velocity v y — vo(x — 0.5) /ro. The taper is defined 
within the band ro < r < n: mass density p = 1 + 9/; x-velocity v x = —fvo(y — 0.5)/r; ^/-velocity v y = fvo(x — 0.5) /r. 

(c) Magnetic field quantities. Across the entire grid: B x = 5; B y = B z = 0. 

The mass density, thermal pressure and magnetic pressure at time t = 0.15 are shown in Fig. [6] In each plot, there are thirty equally- 
spaced contours between the highest and lowest values. These results can be compared with Fig. 18 of Toth (2000) or Fig. 25 of Stone et al. 
(2008), and imply that the magnetic field is divergence free. Again, it is worth noting the 7r-rotational symmetry of the results. 



7 CONCLUSION 

We have presented a partially dimensionally- split algorithm for numerically solving the equations of MHD. This was done by reformulating 
the algorithm of Balsara & Spicer ( 1999 ) and Toth (2000). It is based on the standard reconstruct- sol ve-average strategy (using a Riemann 
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solver), and relies on constrained transport to ensure the magnetic field remains divergence free (V-B = 0). The dimensionally- split approach, 
while not necessarily as accurate as the detailed unsplit algorithm described by Stone et al. (2008), comes with many implementational 
advantages over their approach, and even the Balsara-Spicer-Toth algorithm: the information that needs to be stored during the directional 
passes is minimised; the processing required between complete sets of three directional passes is minimised; the need for staggered grids 
has been completely eliminated (simplifying the use of adapative and non-uniform grids); the algorithm is more easily parallelisable. All 
of these advantages also make the addition of this algorithm to existing dimensionally- split codes much simpler. This makes it particularly 
useful for mature astrophysical codes, which often model more complicated physical effects on top of an underlying hydrodynamics engine 
and so cannot be restructured easily. We included a complete description of the implementation, and illustrative source code will be made 
freely available soon. Finally, we demonstrated the accuracy of the algorithm with several common MHD test problems. 

ACKNOWLEDGMENTS 

This research was supported by the Engineering and Physical Sciences Research Council, and has made use of the resources provided by the 
Edinburgh Compute and Data Facility (ECDF). We would also like to thank Sam Falle at the University of Leeds for his comments. 

REFERENCES 

BalsaraD. S., 1998, Astrophys. J. Suppl. S., 116, 133 

Balsara D. S., Spicer D. S., 1999, J. Comput. Phys., 149, 270 

Colella P., Woodward P. R., 1984, J. Comput. Phys., 54, 174 

Dai W., Woodward P. R., 1998, J. Comput. Phys., 142, 331 

Evans C. R., Hawley J. E, 1988, Astrophys. J., 332, 659 

LeVeque R. J., 1992, Numerical Methods for Conservation Laws. Birkhauser 

Liska R., WendroffB., 2003, SIAM. J. Sci. Comput., 25, 995 

Londrillo P., del Zanna L., 2000, Astrophys. J., 530 

Mignone A., et al, 2007, Astrophys. J. Suppl. S., 170, 228 

Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315 

Orszag S. A., Tang C.-M., 1979, J. Fluid. Mech., 90, 129 

Stone J. M., Gardiner T. A., Hawley J. E, Simon J. . B., 2008, Astrophys. J. Suppl. S., 178, 137 

Toth G., 2000, J. Comput. Phys., 161, 605 



© 2013 RAS, MNRAS 000,[l]{T2 



